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ABSTRACT 


A comparison of a simple recursive ray acoustics algorithm versus a ray acoustics 
algorithm based on solving a system of first-order ordinary differential equations was 
conducted. The recursive ray acoustics (RRA) algorithm was found to be accurate and 
relatively fast. The RRA algorithm is capable of handling sound speed as a function of 
all three spatial coordinates, and this capability was demonstrated. 

Two separate methods of representing a sound-speed profile (SSP) based on data 
points were examined: Akima cubic spline and spatial Fourier series (SFS). The SFS 
representation encountered difficulties in accurately modeling SSPs. Various techniques 
Were applied to improve the SFS sound-speed representation. While accurate sound- 
speed fits were eventually achieved, difficulties remained in the SFS modeling of first and 
second-order derivatives of the sound-speed data. 

The RRA algorithm was tested using the SFS sound-speed representation and found 
to be significantly inaccurate. 

A demonstration was conducted of the ability of the SFS sound-speed represen- 


tation to incorporate randomness in the SSP. 
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I. INTRODUCTION 


The primary focus of this thesis was to test the accuracy and capabilities of a re- 
cursive ray acoustics (RRA) algorithm that can be used to calculate the position, angles 
of propagation, travel time, and path length along a ray path and to generate ray trace 
plots for speeds of sound that are functions of all three spatial coordinates. The RRA 
algorithm is an extension and slight modification of an algorithm briefly discussed by 
Klein [Ref. 1: pp. 29-31]. The algorithm discussed by Klein [Ref. 1: pp. 29-31 ] was ex- 
pressed in terms of an index of refraction while herein, the RRA algorithm is expressed 
directly in terms of the speed of sound. No computer simulation results were presented 
meeiein (Ref. 1: pp. 29-31 |. 

The accuracy of the RRA algorithm was tested by comparing its computer simu- 
lation results with those obtained from a previously developed ray acoustics algorithm 
based on solving a system of first-order ordinary differential equations (ODE) [Ref. 2]. 
The ODE algorithm as applied to ray acoustics was discussed in detail by Lim [Ref. 3]. 
Both the ODE algorithm and the RRA algorithm used the IMSL Version 10 double 
precision computer programm DCSAKM to obtain an Akima cubic spline fit to one- 
dimensional, depth-dependent, sound-speed data for testing and comparison purposes. 
The computer simulation results are presented in Chapter III, Section A. An Akima 
cubic spline was chosen since it is designed to produce a curve that matches the shape 
of the data while minimizing oscillations. In addition, the RRA algorithm was tested 
and shown to be capable of handling a sound-speed profile (SSP) that is a function of 
all three spatial coordinates. This capability was demonstrated by first separating the 
ocean into zones along the z axis. In each zone, the SSP is a different function, in gen- 
eral, of both cross-range x and depth y. A SSP with a gradient in the cross-range di- 
rection will cause a sound ray to leave its initial plane of propagation. Chapter III, 
Section B, includes demonstrations of the special three-dimensional capabilities of the 
RRA algorithm. 

A secondary objective of this thesis was to examine the use of a spatial Fourter series 
(SFS) to fit sound-speed data and to assess its utility relative to an Akima cubic spline 
fit. The appeal of the SFS lay in its representation of spectral content in terms of spatial 
frequencies, as well as its ability to represent derivatives. The precision of the two 


sound-speed fits, as well as that of the first and second-order derivatives of the sound- 


speed data were compared since the first-order derivative as well as the sound-speed fit 
itself are utilized in the ray trace algorithms mentioned above. Chapter IV, Settion A, 
contains the results of this comparison as generated by computer simulation. The SFS 
representation of a SSP was then employed within the RRA algorithm and tested by 
using the same SSPs that were used to compare the ODE and RRA algorithms. Chapter 
IV, Section B, discusses the accuracy of the RRA algorithm incorporating the SFS rep- 
resentacion Of Cie vanous rosie: 

Additionally, random fluctuations in the SSP were modeled using the inherent suit- 
ability of a SFS representation for the inclusion of components at particular harmonics. 
This random element was intended to simulate fluctuations generated by such processes 


as internal waves or by inaccurate collection of sound-speed data. 


If, THEORETICAL ANALYSIS 


A. THE RRA ALGORITHM 
Ifa sound source 1s located at ry = (%, Yo, Z)) and if a ray is launched with initial an- 
gles of propagation f, and @, (see Figure 1 ), then an estimate of the position of subse- 


quent points along the ray path 1s given by (see Figure 2 ) 
ees aera a i= Ls Sanr (2.1) 
where 


Ar, = As; n, (2.2) 


As, 1S an increment of arc length, the magnitude of which is at the user’s discretion, 
ny =u x+vyy tw 2 (2,3) 


is the unit vector along the ray path where 


u; = sin B; cos @,, (2.4) 
v= cos p,, (ee) 

and 
W; = sin B; sin oF (26) 


are the dimensionless direction cosines with respect to the x, y and z axes, respectively, 
and £, and @, are the angles of propagation. 

Since, in general, n,_, is not identical with n,, from Klein [Ref. 1: pp. 29-31] we see 
that, as amplified and applied to the ocean acoustics environment in Ziomek [Ref. 4: pp. 
225 and 230}, 


d A 
oD one, 2.7 





Figure 1. Initial Angles of Propagation f, and @, . 





Figure 2. Ray Geometry. 


Or, approximately, 





A(n(x,)2;) = Vn(x)As; (2.8) 
where n(r,) = 7,, 
n(x) = as (2.9) 


is the dimensionless index of refraction, c(y,) 1s the three-dimensional speed of sound at 
een yy een) (2.10) 
and, c(y,) 1s the three-dimensional speed of sound at 
A AN A 
i ee VAY zy (2.11) 


which 1s the position vector to a point on the ray path after the ith iteration. 


It follows from Equation (2.8) that 





n(r) my — nj) = Vn(r;) As; fl) 
and, as a result, 

A n(r)_ m1 As; ‘ 

a n(r;)) n(y;) Ee vo) 


where, for better accuracy in the evaluation of the gradient, 


NS ie, 


Fat =f) 7 1 (2.14) 


which is the three-dimensional spatial location half-way along the arc length between 
oat and Liy and 


Pe A ft ee: (2.15) 


is the gradient expressed in rectangular coordinates. It should be noted that although 


As, has been shown as a function of iteration number, it will be constant, in general, 
except when it needs to be adjusted to properly handle surface and bottom reflections 
and maximum down range limits. 
In terms of the direction cosines u,v, and w, 
NY )Ui As n(r;_1) 


ya a (275) 


MG) As nr.) 
ce n(x) u n(r;) Oy = 


~ 


and 


— nr, 1) Wy As én(ei_1) 


n(r;) n(r;) dz 





(2.18) 


If c(v,) = c(x,, y, z,) 1s the three-dimensional speed of sound, and if cj—sc(7))miseuae 


speed of sound at the source location, then it follows that 











mee) oe) c(r;) 
ny) (42) Ce (2.19) 
and 
| ente) | dc(r) 
n(r;) Ox =~ «| ae ee hes, (2.20) 


Analogous reasoning applied to the y and z coordinates yields results similar to 


Equation (2.20) , and substitution of these results into Equation (2.13) yields 


c(r;) 


i. C(x) 





A 
Hy 


ni cep] zi veto |, is: (2.219 


In addition, although Klein did not discuss travel time t and path length s calcu- 


lations, these quantities were computed using the following formulae: 


> AS (2.22) 
a es ; 
c(r;_1) 





and 


s=) As. (2.23) 
i= 


Equations (2.1) through (2.23) constitute what 1s referred to herein as the recursive 
ray acoustics (RRA) algorithm first mentioned in Chapter I. Input to this algorithm in 
the form of mathematical formulae for the sound-speed profile or some particular nu- 


merical fit to sound-speed data must be provided to enable the algorithm to function. 


B. SPATIAL FOURIER SERIES REPRESENTATION OF SOUND-SPEED DATA 
Given an ocean D meters deep with a one-dimensional, depth-dependent sound- 
speed profile (SSP), the speed of sound c(y) can be represented by the following spatial 


Fourier series (SFS): 


S 
oy) = leg] +2) ley] cos(2nf,, y — Le), 0<y<D (2.24) 


cS 


where |c,| and 4c, are the magnitude and phase(angle), respectively, of the complex 
Fourier series coefficient c, at harmonic s, |c,| is the magnitude of c, at harmonic s = 0 


and is the “DC component” or average value of the sound-speed data, 


=| 
Sr aa N is 





(2.25) 


is the fundamental spatial frequency in the y direction in cycles per meter where JN is the 
total number of sound-speed data points taken and Ts is the uniform sampling period 
in the y direction in meters per sample, and S is the highest harmonic. Derivatives of 


this Fourier series representation can be directly obtained. 


The complex Fourier series coefficients c, are obtained by taking the spatial discrete 
Fourier transform (DFT) of the sound-speed data as follows: 


V—1 
Came) ens IVE, $= Syne QyoeS (2.26) 
"n= 
where 
HY. = exp( + j2z/.N) (Qa) 


and, when sampling at the Nvquist rate, 
N= 2S + 1. (2.28) 


An attempt to represent a function by a truncated Fourier series will always result 
in oscillations about the functions’ true value. These oscillations normally have small 
enough amplitudes that they are of little consequence. However, an attempt to use a 
SFS representation of a linear SSP with a single constant gradient results in significant 
discrepancies in amplitude near the ocean surface and bottom, due to the inherent dis- 
continuity in the SSP that exists at these boundaries. This difficulty, known as “Gibbs’ 
phenomenon”, does not diminish as more harmonics, that is, sound-speed data points, 
are added to the Fourier series. 

C. Lanczos [Ref. 5] suggested that this phenomenon would be ameliorated if the 
Fourier series could be forced to converge more rapidly by modifying the Fourier series 
coefficients c, by a correction factor, such that the Lanczos smoothed Fourier series co- 
efficients c,, are given by 

C., = — C 70 (2.29) 
where c, is given by Equation (2.25). Hamming [Ref. 6] states that S in Equation (2.29) 
is usually replaced with S + 1 and, thus, our representation of the correction factor will 
be 
sin(s7z/(S + 1)) 


= .30 
Cs, sa(S + 1) C;,, S#O (2.30) 


C. INCORPORATION OF RANDOMNESS INTO THE SOUND-SPEED 
PROFILE 
The speed of sound as a function of depth, c(y) , can be represented as the sum of 


deterministic and random components as follows: 
cly) = cply) + cr), (2.31) 


with cp(y) being the deterministic and c,(y) the random speeds of sound as functions of 
depth y. 


The autocorrelation function of c,(y) is , for any particular y and y’ , 


Re (i y') = Elervery)} . (2.32) 
mere c,(y ) is the complex conjugate of c, at depth y’, O<y'<D. 
it 
N 
cal) = > cy exp( — j2nnfy, y), (2.33) 
n=—N 


for some set of Fourier series coefficients c,, and 
f,, = UD (2.34) 
for D the ocean depth in meters, then 


lies | 
C= j cy) exp( + j2anf, y)dy. (2.35) 


If ca(y) is a random function of depth, then c, will be a random variable. 


Since 
* l DD * 
Elen =e | | ElealecaW’)} exp + rfl —n'v"drdy”, (2.36) 
D* Yo Yo 
where 
E\cely)caly’)} = Rel) (2.37) 


and (c,)° 1s the complex conjugate of c,, and if it is assumed that c,(y) is a wide-sense 


Stationary ( 1.e., statistically homogeneous ) random process, then 
Ro Ay; )') = K, (Ay), (2.38) 
where Ay=y—y’'. If R,, is represented by a finite Fourier series 
N 


Re(Av) =D) tm exP(— J2nmf, Ay), (2.39) 


man-N 


then, as developed in Thomas [Ref. 7: pp. 148-153], 
1 
m= INT Scalttfy,)s (2.40) 


Where S,, 1s the power spectral density of c,(y). 
Thomas [Ref: 7: pp 148-153] goes on to show that 


lara as ion (2.41) 
where d,,, 1s the Kronecker delta function. Therefore, 
2 l 
E} ie, } == are Sealy): (2.42) 
Hence, if c, 1s a zero mean Gaussian random variable with variance r, , that 1s, if 
Cc, = N(0, r,), (2.43) 
then c,(y) will be a zero mean Gaussian random process. 
Tyereiore. 
Reg( sd") = Reg AV) = > tm exP( — j2refy Ay), (2.44) 


n 


which represents the Inverse Discrete Fourier Transform of the power spectral density. 
Thus, knowledge of the power spectral density of c,(y) enables calculation of r, via 
Equation (2.42), which in turn yields c, via (2.43) and cp(y) via (2.33). 
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Il. COMPUTER SIMULATION RESULTS 


A. RAY ACOUSTICS - COMPARATIVE ANALYSIS OF ODE VERSUS RRA 
ALGORITHMS 

Four distinct sound-speed profiles (SSPs) were examined for purposes of compar- 
ative analysis. For aljl cases, a sound source spatial location of x,=0 m and z,=0 m 
was used. The source depth y, was established at a point half-way between the ocean 
surface and bottom, 100 m in most cases. A constant ocean depth of 200 m was mod- 
eled, unless otherwise noted, and for all cases three separate rays were examined. All 
rays were propagated initially in the direction of the positive z axis, that is, 6, = 90°, but 
B, was altered in order to generate unique rays. For SSPs 1 to 3, the launch angle £, 
was set to 45°, 85°, and 135° . SSP 4 was designed to simulate a SOFAR channel. 
Launch angles were therefore chosen to result in the first two rays propagating down 
range without reflection, 1e., 6, = 85° and 95°, while the third launch angle was set to 
By = 135° to test reflection as well. 

The four SSPs considered are defined as follows: 

e SSP 1: Constant speed of sound ( c= 1500 m/s ) 


e SSP 2: Linear SSP with a positive gradient ( c(jy) = c(0) + gy m/s where the speed 
of sound at the surface c(0) = c(y=0) = 1500 m/s and g = +0.017 I/s ) 


e SSP 3: Linear SSP with a negative gradient (c(y) = c(O) + gy m/s where the 
speed of sound at the surface c(0) = c(jy =0) = 1500 m/s and g = —0.017 I/s ) 


e SSP 4: Parabolic SSP with vertex at y, where c(y,) = 1490 m/s and the speed of 
sound at both the ocean surface and bottom is equal to 1500 m/s 


(of = c(y) + aly — Jy)? m/s with a = 0.001 —& for an ocean depth of 200 
im ). 

The above four SSPs generated nine sound-speed data points which were input to 
the Akima cubic spline routine. The results from the ODE algorithm were used 
throughout as the benchmark for comparison purposes. The ODE algorithm used the 
IMSL Version 10 double precision computer program DIVPBS in order to solve the 
system of first-order ordinary differential equations. A range step size of 2.0 m and an 
eunontolerance of | x 10-’ was used in DIVPBS, except for SSP 3, Ray 1! in Table 3, 
where the tolerance had to be reduced to | x 10° in order for the ODE algorithm to run. 
The arc length step size used in the RRA algorithm was 2.0 m unless otherwise noted. 


The results of the various simulation runs are presented in the following seven tables. 
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Table 1. RAY ACOUSTICS RESULTS FOR SSP t, RAYS 1 THROUGH 3 


100.0 | 450 | 100 | 99.95 45.000 9.428091 14.1421 


100.0 oe re 174.89 § 5.000 6.692132 10.0382 


100.0 | 1380 | 100 | 100.00 } 135.000 | 9.428090 | 14.1421 


For SSP 1 (constant speed of sound), Table | shows that at a horizontal range 












r= 10 km, the differences between the two algorithms are as follows: 
e between 0.01 and 0.05 min depth y and 


e 1 usec in travel time t. 


Table 2. RAY ACOUSTICS RESULTS FOR SSP 2, RAYS 1 THROUGH 3 


100.0 162.57 85.487 6.683496 10.0374 


100.0 100.01 | 135.000 | 9.417413 | 14.1421 
| RRA | 1000 135.0 | 100 100.05 | 135.000 | 9.417422 | 14.1421 


For SSP 2 (linear SSP with a positive gradient), Table 2 shows that at a horizontal 













range r = 10 km, the differences between the two algorithms are as follows: 
¢ between 0.04 and 0.07 m in depth y and 


e between 4 and 9 usec in travel time t. 


Table 3. RAY ACOUSTICS RESULTS FOR SSP 3, RAYS 1 THROUGH 3 


100.0 135.0 | 10.0 100.05 | 135.000 | 9.438794 


For SSP 3 (linear SSP with a negative gradient), Table 3 shows that at a horizontal 












range r= 10 km, the differences between the two algorithms are as follows: 
e between 0.04 and 0.07 m in depth y , 
e 0.001 deg in angle of arrival f , and 


® betiveen 5 and 11g seeumetma clecmmacr 


The data presented in Tables | through 3 demonstrate that the two algorithms are 
in very good agreement. Note that the accuracy of the RRA algorithm can always be 


improved by reducing the arc length step size. This is demonstrated in Tables 5 and 6. 


Table 4. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3 


100.0 | sso | 100 37.06 87.244 6.710387 10.0183 


100.0 | 950 | 100 162.94 | 92.756 Ow 10387) 100183 
100.0 Oso eee | Owmuso). | 10.0184 


| ODE 100.0 135.0 | 10.0 144.97 | 134.922 | 9.448852 14.1104 
Paes Pe | ose | ww [eas [oinw [savoir 


For SSP 4 (parabolic SSP), Table 4 shows that at a horizontal range r = 10 km, with 














the ocean bottoin at 200 m, the differences between the two algorithms are as follows: 
e between 0.43 and 0.73 min depthy , 
e between 0.008 and 0.02 deg in angle of arrival f , 
e between 52 and 386 u sec in travel time t, and 


e between 0.1 and 0.6 m in path length s. 


As 1s indicated in Table 4, the differences in travel time are relatively large when 
compared with previous cases. This lower level of accuracy was considered unacceptable 
since, in general, accuracy in the millisecond region 1s normally required for travel time. 
An attempt to improve accuracy was made by reducing the arc length increment used 
in the RRA algorithm from 2.0 to 1.0 m (see Tables 5 and 6). 

The parabolic SSP which produced Table 4 has the most extreme gradient of all the 
SSPs examined and, thus, the two rays launched at 85° and 95° proceed down range 
without reflection from the surface or bottom, as illustrated in Figures 3 and 4. This 
extreme gradient was considered to have been the probable cause of the discrepancies 
noted above and, therefore, the gradient was reduced and the ocean depth increased to 


permit an examination of a similar profile with a more typical gradient (see Table 7). 
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Table 5. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3, 2.0 
M ARC LENGTH STEP SIZE 


100.0 135.0 1230 45.019 4.72441 1 (Ass 


The data in Table 5 for a horizontal range of 5 km and using an arc length step size 













of 2.0 m was generated in order to enable a direct comparison with the arc length step 
size in the RRA algorithm set to 1.0 m, as shown in Table 6. The reduction in final 
range from 10 to 5 km was necessitated by the structure of the RRA computer code, 
which required large quantities of computer storage space to be allocated, particularly 
as arc length step size was decreased. This drawback could be rectified by straightfor- 
ward modifications to the simulation program. 
For SSP 4 (parabolic SSP), Table 5 shows that at a horizontal range of r= 5 km, 

the differences between the two algorithms are as follows: 

¢ between 0.12 and 0.16 min depth y , 

¢ between 0.002 and 0.016 deg in angle of arrival B , 

e between 22 and 95 yu sec in travel time t , and 


e 0.1 min path length s 


Table 6. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS ! THROUGH 3, 1.0 
M ARC LENGTH STEP SIZE 


Foor [wwe | oso | so | oss [sew fasaro [om 
pee [wr |e [or feo | [so 
pe [oa [Ser [eee [os 


For SSP 4 (parabolic SSP), Table 6 shows that at a horizontal range of r=5 km, 









and with the RRA arc length step size set to 1.0 m, the differences between the two al- 
gorithms are as follows: 

¢ 0.06 min depth y, 

e between 0.001 and 0.008 deg in angle of arrival £ , 

® between 11 and 47 yp sec in travel time t , and 


¢ 0.1 min path length s. 


A significant improvement in the consistency of travel time calculations can be seen 
by comparing Tables 5 and 6. Cutting the arc length increment in half reduced the travel 
time discrepancy by approximately one-half. Further reductions in arc length increment, 


an input parameter to the simulation program, would yield further improvements in 
travel time accuracy. 


Table 7. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS t THROUGH 3, 2.0 
M ARC LENGTH STEP SIZE, BOTTOM DEPTH 1000 M 


[o5e | son | sso | wo | rat | vas | araoe [oo 
ran [soo | 50 | woo [as [nas | eras [rom 
Toor [swe | oso | wo | mae | ooo | arse [roms 
rex [soo | 950 1oo [ens | woo | eras [ron 
a 
rex [sn [50 | 100 [ose [ist van ut 


For SSP 4 (parabolic SSP) with vertex c(500) = 1490 m/s, speed of sound at the 
ocean surface c(0) = 1500 m/s, and speed of sound at the ocean bottom c(1000) = 1500 
m/s, i.e., an ocean depth of 1000 m; with c{y) = c(500)+a(y — 500)? where 
a = 0.00004, at a horizontal range r= 10 km, the differences between the two algo- 













rithmis are as follows: 
® between 0.01 and 0.07 m in depth y, 
¢ 0.001 deg in angle of arrival 6 and, 


® between 5 and [5 p sec in travel time t. 


Table 7, when compared with Table 4, both pertaining to an arc length step size of 
2.0 m, demonstrates that a more realistic gradient of 0.02 I/sec enables the RRA algo- 
rithm to handle propagation in a SOFAR channel to a range of 10 km with a much 
improved degree of accuracy. Of course, this could be improved still further by reducing 


the arc length step size. 
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B. RRA THREE-DIMENSIONAL CAPABILITIES 
Figure 5 demonstrates the ability of the RRA algorithm to handle launch angles 


other than @, = 90° . This figure was produced using the parabolic SSP 4 with 


2 


cly) = co) + aly —y), OS ys 200m, (en) 


where 3, = 100m, c(y,) = 1490.0 m/s, a = 0.001 — , and @,—ss . Note that the 
ray stays in the plane created by this angle, as illustrated by the projection of the ray 
onto the Y Z plane as a straight line. 

The RRA algorithm also has the capability to divide the ocean into separate zones 
as a function of z. Each zone can have its own SSP as a function of both cross range x 
and depth y. For computer simulation purposes, the speed of sound was defined as fol- 


lows: 


c(x, y) = clo) + aly SG er = = 200, (3.2) 


| 


where yp = 100m, c(y) = 1490.0 m/s, a= 0.001 => 


and 


OO See = Onin 
&x = (ee) 


=e Ss, — 2s IO Km 


The particular gradient 1n the x direction was selected solely for computer simulation 
purposes to illustrate the capability of the RRA algorithm to accommodate sound speed 
as a function of both x and y. Figure 6 illustrates the ray trace provided by the RRA 
algorithm for the SSP given by Equations 3.2 and 3.3, showing the ray staying in the 
Y Z plane until reaching 5 km, where the gradient in the x direction comes into effect, 


resulting in out-of-plane propagation. 
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IV. REPRESENTATION OF SOUND-SPEED DATA 


A. COMPARATIVE ANALYSIS - SPATIAL FOURIER SERIES VERSUS AKIMA 
CUBIC SPLINE REPRESENTATIONS OF SOUND-SPEED DATA 
Five different exact mathematical functions were used to generate 11 sound-speed 
data points evenly spaced in depth in a 200 m deep ocean. Akima cubic spline and 
spatial Fourier series representations of the SSPs were produced from these data points. 
Percent error calculations of their sound-speed fits as well as their first and second-order 
derivative fits were performed with respect to the Known exact values. If the error be- 
tween the exact and modeled values was less than 10-* , percent error was set to zero. 
The exact functions examined were as follows: 
¢ Constant speed of sound (c = 1500 1n/s ) 
e Linear SSP with a positive gradient ( cy) = 1500.0 + 0.017 y m/s ) 
e Linear SSP with a negative gradient ( cy’) = 1500.0 — 0.017 y m/s ) 


ty 


° ilalf-period sine wave SSP ( c(v) = 1500.0 — 25.0 sin( aa 
¢ Parabolic SSP ( cfy) = 1450.0 + 0.005(v — 100.0)? m/s ) 





) imse 


Tables 8 and 9 with Figures 7, 8 and 9 illustrate the difficulties that the Fourier 
representation encounters in dealing with data derived from an exact constant speed of 
sound. No noticeable difference in the sound-speed ,.it itself is observed for eithiersmne 
Akima cubic spline or the Fourier series representation, but the latter generates some 
small errors in both the first and second-order derivatives. 

Tables 10 and 11 along with Figures 10, 11 and 12 compare the Akima cubic spline 
and Fourier series representations of sound-speed data generated from a linear SSP with 
a positive gradient. In Figure 10 and Table 11 it can be seen that the Fourier represen- , 
tation does not precisely match the sound-speed function, although Figure 10 and Table 
10 show that the Akima cubic spline does. Problems similar to those noted for a con- 
stant speed of sound can be observed by examining Figures I1 and 12 and Table 11. 
Discrepancies in the derivatives are now much larger. 

Tables 12 and 13 along with Figures 13, 14 and 15 show analogous results for a 
linear SSP with a negative gradient. Once again, neither the sound-speed fit nor the 


derivatives are precisely modeled by the Fourier series representation. 
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Constant Speed of Sound. 


Sound Speed Fit - 


Figure 7. 
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First-Order Derivative Fit - Constant Speed of Sound. 


Figure 8 
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Second-Order Derivative Fit - Constant Speed of Sound. 


Figure 9, 
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Sound Speed Fit - Linear SSP with a Positive Gradient 


Figure 10. 
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First-Order Derivative Fit - Linear SSP with a Positive Gradient. 
as 


Figure 11. 
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Figure 12. 


Second Order Derivative Fit - Linear SSP with a Positive Gradient. 
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AKIMA CUBIC SPLINE RESULTS - LINEAR 
SSP WITH A NEGATIVE GRADIENT 
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FOURIER SERIES RESULTS - LINEAR SSP 


WITH A NEGATIVE GRADIENT 
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First-Order Derivative Fit - Linear SSP with a Negative Gradient. 


Figure 14. 
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Second-Order Derivative Fit - Linear SSP with a Negative Gradient. 
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Figure 15. 


The above three exact sound-speed functions were used to compare the Fourier se- 
ries representation with the Akima cubic spline since the Fourier fit can be expected to 
be at its Worst since We are trying to form linear functions from trigonometric ones. The 
next two functions examined, however, are non-linear; the half-period sine wave and the 
parabola. 

Tables [4 and [5 along with Figures 16, 17 and 18 compare the Akima cubic spline 
and Fourier series representations of sound-speed data generated from a half-period sine 
wave. Both techniques encounter some difficulty here, with the Fourier representation 
generally behaving better near the midpoint of the data than the Akima cubic spline fit 
for sound speed and both derivatives. However, near the surface and bottom, the errors 
in the Fourier representation become much larger than those of the Akima cubic spline 
fit. | 

Tables 16 and 17 along with Figures 19, 20 and 21 compare the Akima cubic spline 
and Fourier series representations of sound-speed data generated from a parabola. 
Surprisingly, the Akima cubic spline representation matches sound speed and derivatives 
very well. The Fourier representation behaves in a manner sunilar to its handling of the 
half- period sine wave. It produces small errors near the midpoint of the data and large 
discrepancies at either end for sound speed as well as in the first and second-order de- 
rivauives. 

The disappointing performance of the Fourier series representation of sound speed 
and of the first and second-order derivatives of sound speed led to various attempts to 
salvage it. First, Lanczos smoothing as discussed in Chapter II, Section C, was imple- 
mented. This ameliorated the oscillations of the Fourier representation as shown 1 
Tables 18 and 19 and Figures 22 through 27 for the half-period sine wave and parabolic 
SSPs. It can be seen that Lanczos smoothing results in a smoother fit for sound speed 
and for the derivatives to a point closer to the end points of the data but, nevertheless, 


leaves significant discrepancies near the end points. 
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FOURIER SERIES RESULTS - HALF-PERIOD 


SINE WAVE SSP 
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Figure 16. 
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Sound Speed Fit - Parabolic SSP. 


Figure 19. 
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First-Order Derivative Fit - Parabolic SSP. 


Figure 20. 
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Second-Order Derivative Fit - Parabolic SSP. 


Figure 21. 
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Wave SSP. 


Figure 24. 





@ectewmewewewe ee weweweeweaxe ae & 


SNIWdS JISND YWIAY 
TH00W O33dS-ONNdS 
SIN{Od UIUO @ 

ON3971 








ONIHIOOKS SGZONUT HLIM 
OfYSLIU ONY 3Y0s3@ SINIOd ULUG OILSHINAS 40 YSGWIN TEISINIOd ULYO JO Y3GHN TENISINO 


dSS O170GuuUd 


i 
: 
' © cweccevacscerccete > ocewceweccecesccece oat 
i 
2 ; 
; : 
[tecnesscereornees t pena sean eaeeeeeae oat 
| i 
Seeecceecwseeeroeseot See Seeeees cs eeettece Ovt 
sendsaceaeiesannaes ee eee oe eeeneeeennaee Ort 
: 
; 
: 

*feeee eeacessoees o eeetcoe oot 


coeveo 
ae 


‘ deeun teseueseenaece 4 eaeseveneeeenes aves oa 
: 
[eveeressenctecnnnne Focceeeeereseseneees os 
i : 
: : 

eeeeeeee eee steeceet coneeccccescccnnees mecsccteccccnconecce oP 
i 
i eee none eeaen + cane penesseneersoes rm 
{ ; 
. Q 





(AJO JO NOILY OdYSINI WOLYSWAN YO3 SSNOINHOSL 


Sound Speed Fit with Lanczos Smoothing - Parabolic SS P. 


Figure 25. 
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First-Order Derivative Fit with Lanczos Smoothing - Parabolic SSP. 


Figure 26, 
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Figure 27. 
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The difficulty in achieving a good fit at or near the end points of the data using a 
spatial Fourier series representation is a result of the instantaneous commencement and 
cessation of sound-speed data at the ocean surface and bottom. Therefore, an attempt 
to overcome this difficulty was made by adding synthetic data points above the ocean 
surface and below the ocean bottom so as to enable the Fourier series to “settle down’ 
before the real data begins. 

Synthetic data was created by first examining the three data points nearest the ocean 
surface and bottom. Similar procedures were used to generate synthetic data both above 
the surface and below the bottom. Therefore, only the surface case will be examined. 

lect 


c(v,) — c(0) 
(a (4.1) 
where c(0) = c(y=0) , and 
c(y) — cQ) 


where y, is depth in meters at the ith data point. If we let 
Am = th, — Myo (4.3) 
and if 
|Arz| = |p, — mol <eé (4.4) 
for some arbitrarily small e, then linear extrapolation is used, L.e., 
Gee — CO) + mij y_;, (= 1,2,..,9YNPTS (4.5) 
where 
=| ne (4.6) 


and SYNPTS is the number of synthetic data points desired. 


ao 


If the test at Equation 4.4 fails, and if 
| tig... << aleeten le (4.7) 


then mirror imaging about the ocean surface is appropriate, L.e., 
c(y_,) = eQy), ¢ = 12),5 eee (4.8) 


If both of the tests at Equations 4.4 and 4.7 fail, then neither linear extrapolation 
nor murror imaging 1s appropriate. Under these circumstances it is appropriate to re- 


verse the rate of change of the slope and let 
c(y_,;) = c(O) + (my — (Am)y_,, i = 1,2,..,SYNPTS . (4.9) 


This methodology of extending the real data was implemented using five synthetic 
data points both above the surface and below the bottom. The results are shown in 
Tables 20 through 24 and Figures 28 through 42. 

A substantial improvement in the accuracy of the Fourier series representation can 
be seen with respect to all the functions examined. Nevertheless, some difficulty 1s still 
evident near the ends of the original data. A final attempt to resolve this discrepancy 
was made by using the Akima cubic spline representation of sound speed as input to the 
Fourier algorithm, and sampling the SSP generated by the Akima cubic spline fit at some 
arbitrarily high frequency. This was attempted using a sampling rate ten times greater 
than available from the original data. The synthetic data as described above was not fit 
using the Akima cubic spline method. Results for the half-period sine wave and 
parabolic SSPs are presented in Tables 25 and 26 and Figures 43 to 48. 

These results are once again an improvement on previous work. Nonetheless, some 
inaccuracies remain. Unfortunately, time constraints precluded further, more sophisti- 


cated attempts to further improve the accuracy of the Fourier representation. 
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Sound-Speed Fit - Constant Speed of Sound with Synthetic Data. 


Figure 28. 
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(WIHL&30 
First-Order Derivative Fit - Constant Speed of Sound with Synthetic 


Data. 


Figure 29. 
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Data. 


Figure 30. 
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Figure 31. 
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Figure 36 
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FOURIER SERIES RESULTS WITH LANCZOS SMOOTHING, 
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Sound-Speed Fit - Half-Period Sine Wave SSP, Akima Cubic Spline 


Input to Fourier Method. 


Figure 43. 
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First-Order Derivative Fit - Half-Period Sine Wave SSP, Akima Cubic 


Spline Input to Fourter Method. 


Figure 44. 
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Second-Order Derivative Fit - Half-Period Sine Wave SSP, Akima Cu- 


bic Spline Input to Fourier Method. 


Figure 45. 
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Figure 46. 
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Derivative Fit - Parabolic SSP, Akima Cubic Spline Input 


First-Order 
to Fourier Method. 


Figure 47, 
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Second-Order Derivative Fit - Parabolic SSP, Akima Cubic Spline Input 


Figure 48, 
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B. RAY TRACE RESULTS USING SPATIAL FOURIER SERIES 
REPRESENTATION OF SOUND-SPEED DATA WITH THE RRA ALGORITHM 
1. Sound Speed Including Deterministic Component Only 

The RRA algorithm using the simple Fourier series sound-speed representation 
including Lanczos smoothing but without synthetic data or Akima cubic spline input to 
the Fourier method was tested using the same source locations, SSPs, and rays as in 
Chapter III, Section A. That is, the SSPs were once again defined as follows: 

e SSP 1: Constant speed of sound (c = 1500n1/s) 


e SSP 2: Linear SSP with a positive gradient (c{y) = c(O) + gy m/s, where the 
speed of sound at the surface c(0) = c(y=0) = 1500 m/s and g = +0.017 1/s ) 


© SSP 3: Linear SSP with a negative gradient (c(y) = c{(O) + gy m/s, where the 
speed of sound at the surface c(0) = ciy=0) = 1500 m/s and g = —0.017 I/s) 


e SSP 4: Parabolic SSP with vertex at y, where c(y,) = 1490 m/s and the speed of 
sound at both the ocean surface and bottom is equal to 1500 m/s 


- = c(y,) + a(y — 100)? m/s with a = 0.001 1 for an ocean depth of 200 
ia) ): 

The sound source spatial location was established at x, =0 mandz,=0Om. A 
source depth of y, = 100 m was established for an ocean depth of 200 m, unless otherwise 
memed. Launch angles 8, were 45°, 85° and*135° for SSPs ly2-and 3; and 85°, 95°, and 
135° for SSP 4, which again was designed to simulate a SOFAR channel. 

The ODE algorithm using the Akima cubic spline fit to sound-speed data re- 
mained as the benchmark for comparison purposes, as detailed in Chapter III, Section 
A. For ease of reference, the results of the ODE algorithm are reproduced in the tables 
of this sub-section. The ODE’s range step size was left at 2.0 m and the arc length step 
size in the RRA algorithm was left at 2.0 m as well, except where otherwise noted. 

The results of the various computer simulation runs are presented in Tables 27 
to 33 below. 
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Table 27. RAY ACOUSTICS RESULTS FOR SSP 1, RAYS 1 THROUGH 35 
FOURIER SOUND-SPEED FIT 


Trex [we | 0 |e | ee | 0 | va | 
reer [on [se _[ or [ves [so [coum [no 
Pe [ oe | | on [| Sn 
Tee rom [ne a Peer 


For SSP | (constant speed of sound), Table 27 shows that at a horizontal range 
of r= 10 km, the differences between the benchmark ODE algorithm and the RRA al- 


gorithm using the Fourier series sound-speed fit are as follows: 











¢ between 0.01 and 0.06 m in depth y and 


e 2 msec in travel time t 


Note the significant degradation of the RRA results relative to those obtained 
using the Akima cubic spline fit as in Table 1 of Chapter HII, Section A. This compar- 
ative degradation of consistency between the ODE and RRA algorithms will be evident 


throughout this subsection. 
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Table 28. RAY ACOUSTICS RESULTS FOR SSP 2, RAYS 1 THROUGH 3, 
FOURIER SOUND-SPEED FIT 


Foor [vo | so | we om [son [saree [mmr 
Trex [oo | 50 | wo [over [som [saan [mr 
Toor [we | aso | wo [rae [var | canoe [wore 
EC 
Fone [me | we | wo frm [oso] var [rarer 
Paes [rma | sso | 0 [ms [ros var [ae 


For SSP 2 (linear SSP with a positive gradient), Table 28 shows that at a hori- 















zontal range r = 10 km, the differences between the benchmark ODE algorithm and the 
RRA algorithm using the Fourier series sound-speed fit are as follows: 
¢ between 0.02 and 0.19 m in depth y, 
e (0.005 deg in angle of arrival £, 
e between 16 and 58 wu sec in travel time t, and 


e 0.1 min path length s. 
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Table 29. RAY ACOUSTICS RESULTS FOR SSP 3, RAYS 1 THROUGH 3, 
FOURIER SOUND-SPEED FIT 


100.0 | 450 | 100 | 99.99 45.000 9.438784 14.1421 


100.0 168.89 | 84.512 | 6.699404 | 10.0380 


For SSP 3 (linear SSP with a negative gradient), Table 29 shows that at a hon- 













zontal range r = 10 km, the difference between the benchmark ODE algorithm and the 
RRA algorithm using the Fourier sound-speed fit are as follows: 
¢ between 0.03 and 0.13 min final depth y, 
¢ (0.020 deg in angle of arrival f, 
e between 14 and 67 u Sec in travel time t, and 


e 0.1 min path length s. 


The results contained in the above tables, while poorer than those detailed in 
Chapter III, Section A, nevertheless produce travel times which are consistent in the 


millisecond region. However, this is not the case for the SSPs examined next. 
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Table 30. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3, 
FOURIER SOUND-SPEED FIT 


Foe [wwe | aso | ww | om | oa | crosr [ mor 
es 
Tres [oo [wo | wo [ne [ns [ons | mons 
Toor Pawo ose | we [mn [oom same [mi 
Pas [me [oso | we [189s [pnw [sass [a 


For SSP 4 (parabolic SSP), Table 30 shows that at a horizontal range of r = 10 
km, the differences between the benchmark ODE algorithm and the RRA algonthm 


using the Fourier series sound-speed fit are as follows: 









¢ between 0.99 and 5.8 m in final depth y, 
e between 0.008 and 0.570 deg in angle of arrival £, 
¢ between 337 and 654 wu sec in travel time t, and 


e 0.2 min path length s. 


As in Chapter III, Section A, SSP 4 shows significantly less consistency be- 
tween the two algorithms than the previous SSPs examined. Again, two approaches 
were taken to attempt to reduce this discrepancy. First, the arc length step size was re- 
duced from 2.0 to 1.0 m (see Tables 12 and 13) and second, the steepness of the gradient 
was reduced from 0.1 to 0.02 1/sec (see Table 14). 
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Table 31. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3, 2.0 
M ARC LENGTH STEP SIZE, FOURIER SOUND-SPEED FIT 


For SSP 4 (parabolic sound speed), Table 31 shows that at a horizontal range 









of r = 5 km, using an arc length step size of 2.0 m, the differences between the 
benchmark ODE algorithm and the RRA algorithm using the Fourier series sound-speed 
fit are as follows: 

¢ between 0.32 and 4.54 m in final depth y», 

e between 0.002 and 0.162 deg in angle of arrival £, 

¢ between 188 and 249 yu sec in travel time t, and 


e 0.2 min path length s. 


The data in Table 31 above was generated in order to enable a direct compar- 
ison with the arc length step size set to 1.0 m, as shown in Table 32 below. The re- 
duction in final range from 10 to 5 km was necessitated by the structure of the RRA 
computer code, which required large quantities of computer storage space to be allo- 


cated, particularly as arc length step size was reduced. 
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Table 32. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3, 1.0 
M ARC LENGTH STEP SIZE 


[oe | we | aso | so [ms [oon | somrm me 
Tas [mo | aso | so | om | oe | sass | som 
fren [oe | aso | s0_[ one [ass [sen [sm 


For SSP 4 (parabolic SSP), Table 32 shows that at a horizontal range of r= 5 
km, and with the RRA arc length step size set to 1.0 m, the differences between the 
benchmark ODE algorithm and the RRA algorithm using the Fourier series sound-speed 
fit are as follows: 
















e between 0.24 and 4.61 m in final depth y, 
e 0.155 deg in angle of arrival £, 
® between 175 and 285 usec in travel time t, and 


e¢ 0.2 min path length s. 
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Table 33. RAY ACOUSTICS RESULTS FOR SSP 4, RAYS 1 THROUGH 3, 2.0 
M ARC LENGTH STEP SIZE, BOTTOM DEPTH 1000 M, FOURIER 
SOUND SPEED FIT 


100.0 768.46 Od heed 6.705420 10.0152 


RRA 100.0 231.54 | 86.478 6.705420 | 10.0152 


100.0 135.0 545.01 | 134.997 | 9.448829 | 14.1104 
100.0 544.58 | 134.997 | 9.448180 | 14.1107 


For SSP 4 (parabolic speed of sound), with vertex c(500) = 1490 m/s and speed 














of sound at the ocean surface c(0) = 1500 m/s and speed of sound at the ocean bottom 
c(1000) = 1500 m/s, ie., an ocean depth of 1000 m; at a horizontal range of r= 10 km, 
the differences between the benchmark ODE algorithm and the RRA algorithm using 
the Fourier series sound-speed fit are as follows: 

e between 0.43 and 8.95 m in final depth y, 

e 0.134 deg in angle of arrival £, 

e between 478 and 649 yu sec in travel time t, and 


e between 0.2 and 0.3 m in path length s. 


Table 33 shows no significant improvement over Table 30 for SSP 4 with the 
original extreme gradient. The errors associated with the Fourier sound speed repre- 
sentation have resulted in this lack of improvement, whereas substantial improvement 
was observed between these two cases for the Akima cubic spline fit with the RRA al- 
gorithm (see Tables 4 and 7 in Chapter III, Section A). 
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2. Sound Speed [Including Deterministic and Random Components 
In order to investigate the effects of randomness in the SSP, a mode] of the 
power spectral density of the random component of the speed of sound c,(y) had to be 
proposed. The most straightforward model of band-limited white noise was used, as 


given by 


Somme ec BEV 
S = a 4.10 
cay) ‘o lee 2 HY em 


where BW is the bandwidth, and S, 1s the constant magnitude of the power spectrum at 
all spatial frequencies f, . 
AS a consequence of the assumption of zero mean white noise, the variance of 


caly) is given by 





on, = 2S, BW (4.11) 
and, as a result, 
oo. 


A value for the variance of c,(y) in shallow water of 2.5 x 10-’ was obtained 
from Brekhovskikh and Lysanov [Ref. 8: p. 208 |. The simulation program was run with 
randomness incorporated as above and for SSP | (constant speed of sound). When the 
results were compared with those from the RRA algorithm previously run using the 
Fourier series representation of sound-speed data, differences of a maximum of 0.1 min 
final depth y, 0.002 deg in angle of arrival f, and 2 yw sec in travel time t were observed. 
The RRA algorithm was then run with SSP 4 (parabolic SSP in a 200 m deep ocean) 
incorporating randomness and differences of a maximum of 0.22 m in final depth y , 
0.018 deg in angle of arrival B, and 16 » sec in travel time t were observed. 

The variance was then increased to determine at what magnitude it caused the 
randomness to become significant. Further simulation runs revealed that for SSP I, 
differences in travel times of greater than 10 4 sec needed o?, to be at least of the order 
of 10- , four orders of magnitude greater than the estimate given in Brekhovskikh and 
Lysanov [Ref. 8: p. 208 ]. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The recursive ray acoustics (RRA) algorithm was shown to be an accurate ray 
acoustics algorithm when used in conjunction with an Akima cubic spline fit to sound- 
speed data. The RRA algorithm can be made as accurate as desired by reducing the arc 
length step size. Although the ODE algorithm computes amplitudes which the RRA 
algorithm does not, the RRA algorithm was observed to be significantly faster in terms 
of computer eExccution cme: 

The abilitv of the RRA algorithm to handle sound speed as a function of all three 
spatial coordinates enables its application to real world problems wherein a horizontal 
gradient of sound speed exists, either in the cross-range or down-range direction. This 
could conceivably be the case for either long range ray propagation or for propagation 
in estuaries Or across current boundaries where significant horizontal gradients can exist. 

A spatial Fourier series (SFS) representation of sound-speed data, while attractive 
for reasons of its spectral content and the ease of incorporation of randomness at spe- 
cific harmonics, encountered difficulty with regard to providing accurate fits to sound- 
speed data. Several techniques were used to attempt to rectify this, including the use 
of Lanczos smoothing , addition of synthetic data above the ocean surface and below 
the ocean bottom, and sampling of the Akima cubic spline fit to provide data points fer 
input to the SFS method. These techniques resulted in some substantial success in the 
reduction of the inaccuracies of the SFS method. However, although the SFS method 
eventually provided accurate fits to the sound-speed data, it still provided inaccurate first 
and second-order derivatives when compared against the results obtained from the 
Akima cubic spline method. 

A comparison was made between the ODE algorithm using an Akima cubic spline 
sound-speed representation and the RRA algorithm using a SFS sound-speed represen- 
tation, with the SFS representation in its most basic form, that is, with none of the 
techniques referred to in the above paragraph incorporated. The RRA/SFS combina- 
tion was found to be wanting in accuracy. 

As was mentioned previously, the RRA algorithm using the Akima cubic spline fit 
to sound-speed data was proven fast and accurate, and can be made as accurate as de- 


sired by reducing the arc length step size. However, investigation is warranted to de- 
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termine if the RRA algorithm can be generalized to include amplitude calculations, 
especially at focal points and, hence, caustics. 

Throughout the simulations undertaken in this thesis, a zero-slope, flat ocean sur- 
face and bottom was assumed. In order to apply the RRA algorithm to realistic ocean 


boundaries, the ability to accommodate arbitrarily shaped surfaces should be included. 
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